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ABSTRACT 


This thesis investigates the use of wavelet transforms in the detection and 
estimation of spread spectmm frequency hopping signals. The technique developed in 
this work makes only two basic assumptions of a minimum hopping time and a minimum 
frequency hopping differential. The approach is based on the phase information of the 
temporal correlation function and the resulting discrete wavelet transform is used to 
estimate the hopping time of frequency hopping signals. Results show the proposed 
scheme is robust to additive white noise for SNR levels of 3 dB and above. 
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1. INTRODUCTION 


A. BACKGROUND 

Spread spectrum communications schemes have received ever increasing attention 
over the past two decades as numerous civilian applications have joined the military 
applications [1]. Originally, spread-spectrum communications were developed to provide 
a secure means of communication in the hostile environments of World War n [1]. This 
work remained mostly classified until 1970s, but by the late 1970s the literature was 
starting to amass [1]. Naturally, it was not long before civilian applications for spread 
spectrum communications began to develop. One such system, the Global Positioning 
System (GPS), which was originally developed for the military, soon found widespread 
civilian application in position location for civilian vehicles, commercial vehicles and 
ships, and hunters and fishermen [1]. Other applications include cellular telephony and 
personal communication systems, both of which are still growing rapidly [1]. 

With this increased use of spread-spectrum communications, naturally, came the 
question of unauthorized interception. This question was not only of importance to the 
military, but also to communications regulatory boards, such as the Federal 
Communications Commission (FCC) [2]. This question has, therefore, given rise to a 
considerable amount of research over the past decade focused on the detection and 
estimation of spread-spectrum communications signals [2]-[6]. Two main assumptions 
typically found in the literature are that the hop timing is constant and known, and that the 
hopping frequencies are selected from a known class of candidate frequencies. Even 
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when the hop timing is not assumed known, it is still usually assumed constant [3]. 
These assumptions generally restrict the detection and estimation schemes to frequency 
hopping (FH), one of the more popular spread-spectrum communications techniques. 

During approximately the same period as the detection and estimation of FH 
signals research was being conducted an analysis technique called wavelet analysis began 
finding wide spread use in signal processing [7]. One reason wavelet analysis was of 
particular interest is its multiresolution analysis capabilities. In other words, wavelet 
analysis allows one to “zoom in” for a detailed look at a signals characteristics, or “zoom 
out” for a global view of the signal. To use a map analogy, the former would be a map of 
city streets (a “large” scale map) while the latter would be a map of the United States (a 
“small” scale map). 

A spread-spectrum communications signal consists of both long duration 
relatively low frequency components (“large” scale) in the carrier frequencies and short 
duration high frequency components (“small” scale) in the transitional hops. As 
mentioned above, wavelet analysis is, in theory, well suited for analyzing signals of this 
type. Indeed, the two have met [4] and wavelet analysis was shown effective in detecting 
frequency hopping signals. 

B. GOALS 

The primary goal of this thesis is to provide a new approach for the detection and 
estimation of frequency hopping signals which makes none of the restrictive assumptions 
listed above. By not making such restrictive assumptions, it is hoped that a secondary 
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goal of wider application to the detection and estimation of other spread-spectrum 
communications techniques (i.e., direct sequencing, time hopping and hybrids of the 
three) given in Peterson [1] can be obtained. 

This thesis is composed of six chapters with this introduction being the first. 
Chapter n introduces frequency hopping signals which are derived from one particular 
spread-spectrum communications signaling technique. Chapter lH defines and explains 
wavelet analysis by analogy with the more familiar Fourier analysis. Chapter IV provides 
the theory behind wavelet-based transient detection and introduces a few preprocessing 
tools which will be used in the detection and estimation algorithm. Chapter V 
enumerates the steps of the algorithm, describes the simulation used in the testing of the 
algorithm, applies the algorithm to an example FH signal, and provides the results of the 
simulation. Chapter VI provides a summary, conclusions and proposed further study. 
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II. FREQUENCY HOPPING SIGNALS 


In communications there exist techniques for spreading the spectrum of 
transmitted signals. These techniques are called spread spectrum (SS) because the actual 
transmission bandwidth is much greater than the minimum bandwidth which would 
otherwise be required to transmit the given information [8]. Such techniques provide 
several benefits which include signal interference suppression, tow probability of signal 
detection, and multiple access by many signals to the same spectral region. This section 
introduces frequency hopping (FH) signals which are one candidate scheme in SS 
applications. We first describe the concept of frequency hopping, then give a simple 
example, and finally, discuss some of the aforementioned benefits of spectrum spreading. 

A. SPREADING THE SPECTRUM 

It has been shown that for signals of bandwidth W and duration T that the 
dimensionality of the signaling space is approximately equal to 2WT [9]. The idea 
behind SS is to increase this dimensionality (i.e., to increase the size of the signal space), 
in order to realize the benefits mentioned above. It can readily be seen that there are only 
two options for accomplishing this spread, either increase W or increase T. One method 
for increasing the bandwidth, W, is called frequency hopping. Increasing the time 
variable. T, can be accomplished via time hopping (TH). A detailed survey of other SS 
techniques and/or TH techniques can be found in Sklar [8]. 
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Figure 2.1: FH/MFSK system 

Although the particular modulation scheme is not of great importance to the 
detection algorithm described in this thesis, M-ary frequency shift keying (MFSK) will be 
used in the discussion and example which follow since this modulation scheme is most 
commonly found in FH systems. It should be noted, however, that it will be the FH 
properties of the signal, not the particular modulation scheme which will be of 
importance later when the detection algorithm is described. In FH systems the available 
bandwidth, W**, is subdivided into a large number of frequency slots, N. The data symbol 
modulates the carrier frequency, fk, where k is selected by a pseudorandom number (PN) 
generator to be between 1 and N. The FH system can be thought of as a two-step 
modulation process, although in practice the two steps would be combined into one. The 
first step. MFSK modulation, would be followed by the second step. FH modulation, as 
shown in Figure 2.1. Note that the pattern of transmitted signal in the time-frequency 
plane will be affected primarily by the FH modulation step, regardless of the particular 
modulation scheme or the number of symbols transmitted per hop. For a given hop, the 
occupied bandwidth, W, is the same as that for conventional MFSK. However, averaged 
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over many hops a much greater bandwidth equal to Wss is occupied, and the spectrum has 
been spread. 

B. FREQUENCY HOPPING EXAMPLE 

Assume that the bandwidth, Wss, is equal to 20 MHz (5 MHz - 25 MHz) and that 
the minimum frequency hopping interval, Af, is equal to 1 KHz. This implies that fi=5 
MHz, f 2 = 5.001 MHz, f 3 = 5.002 MHz, ...,fN=25 MHz where: 

W 

= 77 , ( 2 . 1 ) 

A/ 

equals 20,000. Now, if the PN generator produces the sequence {7500, 1250, Mill, 
11667, 2143}, the first hopping frequency will be: 

/7500 = 5,000,000 + A/ • 7,500 = 1 15 MHz, (2.2) 

and similarly, the next four frequencies will be fi 25 o= 6.25 MHz, fi 7727 = 22.727 MHz, 
fii 667 = 16.667 MHz, and f 2 i 43 = 7.143 MHz. The pattern for this particular FH signal in 
the time-frequency plane is shown in Figure 2.2. For a given hop (for example, between 
0 and 20 |is) we see that the occupied bandwidth, W, is the same as that for conventional 
MFSK. However, averaged over many hops, a bandwidth of W** is occupied. Hence, the 
spectrum has been spread from W to Wj*. Note that depending on the definition of 
bandwidth that this example does not completely cover Wjs. However, it should be 
obvious that as transmission time increases and the carrier frequency, fk, for each hop is 
randomly selected, eventually, the entire W^s will be covered. 
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Figure 2.2: An example of a FH pattern 


C. BENEFITS OF SPREAD SPECTRUM 

1. Interference Suppression 

Consider signals of bandwidth Wss and duration T where Wjs is the SS bandwidth. 
By definition, white Gaussian noise (WGN) has infinite power spread uniformly over all 
frequencies (8]. It can be seen from this definition that regardless of where within W^s 
the actual information resides, it suffers from the same amount of interference from 
WGN. However, jammers must be of a fixed finite power. Therefore, the jammer must 
decide to either spread this power in small quantities over the whole or a large portion of 
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Figure 2.3: Effect of spectrum spreading (a) in the presence of GWN and 
(b) in the presence of an intentional jammer. From Ref. [8] 

Wss or to concentrate this finite power on one or a few smaller portions of the total 
bandwidth, Wss- 

Figure 2.3 shows the effects, in the frequency domain, of spreading the spectrum 
both in the presence of GWN and an intentional jammer. In Figure 2.3a, the shaded 
region. No, is the piower spectral density (PSD) of GWN, G(f) is the power spectrum of a 
signal g(t) before spreading it from signal space W to Wss, and Gss(0 is the power 
spectrum of g(t) after spreading its signal space. We can see that interference due to 
GWN has the same effect on the signal both before and after the spectrum spreading, 
since the GWN’s infinite power is spread uniformly over all frequencies. In Figure 2.3b, 
J is the finite |X)wer of the jammer and the shaded region, J’o = JAV, is the jammer’s 
power spectral density. Here we see, in contrast to the GWN case, that SS does reduce 
the interference due to an intentional jammer, since it typically will not occupy the same 
spectral location as the jammer for the duration of the transmission. Once the spectrum is 
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spread, the jammer has the choice of either covering the entire spread spectrum 
bandwidth, Wss, lightly with its finite power as shown in Figure 2.3b (middle), or of only 
covering only a portion of Wss with increased concentration of its finite power as shown 
in Figure 2.3b (bottom). Note how the shaded region which is the PSD, Jo, of the 
jamming signal is, in the case of Figure 2.3b (middle), wide and short, indicating that the 
effective jamming power is distributed relatively lightly over all of Wss- In the case of 
Figure 2.3b (bottom), however, the shaded region is narrower, but taller, indicating that 
the effective power being delivered is increased by reducing the spread spectrum signal 
bandwidth coverage to, p (0 < p < 1), a portion of Wss. Therefore, the jammer has 
chosen to concentrate its finite power, J, to a portion of the spread spectrum, Wss , in 
order to have an effect on the SS signal when it resides in this portion (the shaded region 
of Figure 2.3b (bottom)). 

2. Energy Density Reduction 

Another benefit of spread spectrum techniques in signaling is that they provide a 
low probability of detection (LPD) [8]. As seen earlier a SS signal’s energy is spread 
throughout a bandwidth much larger than that used in conventional schemes. Therefore, 
the signal energy pre.sent at any given frequency is very small. Hence, the SS signal will 
appear buried in noise to any potential receiver which does not possess the synchronized 
spreading signal for FH demodulation. 

Energy density reduction also is desirable in unrelated and broader (i.e., civilian as 
well as military) applications. In satellite communications, the downlink transmissions 
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must meet international regulations on the spectral density which is transmitted [8]. 
Spread spectrum schemes allow an increase in the total transmitted power while adhering 
to these regulations. 

3. Multiple Access 

A benefit of spread spectrum signaling which has found wide application in 
cellular telephony is code-division multiple access (CDMA) [1]. This technique provides 
simultaneous users with separate unique spreading signal codes. Therefore, many 
simultaneous users can share the same spectrum, Wss, by simply not occupying the same 
portion of the spectmm at the same time. The insurance that they will not occupy the 
same coordinates in Wss is provided by the unique spreading signal code. In addition, 
CDMA offers some measure of privacy as it will be difficult to decode the signal without 
a synchronized replica of the spreading signal code. 

In this chapter we have introduced/re^Mcncv hopping signals. Particular attention 
was paid to the FH properties of the signal instead of the particular modulation scheme, 
since ihc.se arc the properties which will be used to detect and estimate hopping times. In 
the following chapter wavelets and wavelet transforms are introduced. The wavelet 
transform will be the primary tool u.sed in our detection and estimation algorithm 
developed in Chapter V. 
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III. WAVELETS 


Wavelet analysis has become widely known within the last ten to twelve years, 
although similar analysis techniques were employed in various disciplines within 
engineering, physics, and mathematics as early as the beginning of the century [7]. There 
are various methods for describing and explaining wavelet analysis. However, the most 
useful method might be that based on analogies with the more familiar classical Fourier 
analysis [7, 10, 11]. Thus, first, we will provide a brief review of Fourier analysis, and 
next present an introduction to wavelet analysis. 

A. FOURIER ANALYSIS 

1. Fourier Series 

a. Continuous-Time Representation 

Recall that a periodic function, x(t), with a fundamental period To and 
fundamental frequency /o = ^ , has a complex Fourier series expansion in terms of 

complex exponentials given by: 

^(/)= (3.1) 

where e are the set of harmonically related complex exponentials. 
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The Fourier series coefficient, a*, represents the amount the A:-th harmonic of the basis 
functions, e , is present in x(t). Multiplying both sides of Equation 3.1 by e 
and then integrating the result over one period To, leads to: 




(3.2) 


t*T 

Note that dt - Q for k n, as the integration is carried out over one 


full 


period of the complex exponential. Thus: 





Therefore, Equation 3.3 leads to: 


(3.3) 


= dt . (3.4) 

''O T„ 

Equations 3.1 and 3.4 define the Fourier series of a periodic signal. Recall that Equation 
3.1 is called the synthesis equation and Equation 3.4 the analysis equation [12]. Note that 
the ao coefficient given as: 


^0 “ 7 " ' 

^0 T„ 

represents the dc component, or average value, of x(t). 


(3.5) 
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b. Discrete-Time Representation 

For discrete time periodic applications, the discrete time Fourier series 
expansion of a periodic signal, x[n], is given by: 

k=0 

where the Fourier series coefficient, a*, is given by 

n=0 

2. Fourier Transform 

a. Continuous-Time Representation 

The introduction of the Fourier transform which allows the representation 
of not only periodic (as we saw in the last section), but also aperiodic signals as linear 
combinations of complex exponentials was one of Fourier’s most significant 
contributions to this form of signal analysis [12]. 

The Fourier transform is obtained by examining the limiting behavior of 
the Fourier series representation as the period. To, is allowed to grow arbitrarily large to 
infinity [12]. The resulting Fourier transform pair is given by: 


x(r)= 

(3.8) 


(3.9) 


(3.6) 


(3.7) 
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By convention, the analysis equation, Equation 3.9, is called the Fourier transform, while 
the synthesis equation. Equation 3.8, is called the inverse Fourier transform [12]. Note 
that the Fourier series coefficients can be calculated via the Fourier transform since, for 
periodic signals, X(f) will have amplitudes {ak} as given by Equation 3.4 which will occur 
only at a discrete set of harmonically related frequencies, l<f, k=0, ±1, ±2, .... In contrast, 
for aperiodic signals the complex exponentials of Equation 3.9 occur with amplitude 
X(f){df) over a continuum of frequencies [12]. 

b. Discrete-Time Representation 

Fourier synthesis and analysis equations can also be defined for discrete 
time aperiodic signals using a discrete-time representation of the Fourier transform. 
Letting Q = 27tf represent the digital frequency to distinguish from the continuous 
frequency /, leads to: 

x[n] = fx(Q)cda, (3.10) 

X(a)= (3.11) 


3. Short-Time Fourier Transform 

While the Fourier transform works very well for stationary signals, its limitations 
quickly become apparent when this analysis technique is applied to a non-stationary 
signal. Consider the signal from the example in section E-B. The signal Xi(t) consists of 
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one frequency, f 75 oo = 12.5 MHz, for 0 < t < 20 ^is and then a hop to frequency, fi 25 o = 
6.25 MHz, for 20 |Lis < t < 40 |is, as shown in Figure 3.1a. The Fourier transform Xi(/) 
clearly shows the two component frequencies of this signal, as illustrated in Figure 3.1c. 
Now, consider a second signal, X 2 (t), of a duration equal to 40 |xs. However, instead of a 
having a frequency hop the signal is now made up of the summation of two pure 
sinusoidal signals X 75 oo(t) and xi 25 o(t) with frequencies, f 75 oo and fi25o, respectively, each 
with duration of 40 |is, as shown in Figure 3.1b. The resulting Fourier transform X 2 (f) is 
very similar to X\(f), as shown in Figure 3.Id. Thus, this example illustrates the fact that 
the Fourier transform does not provide any temporal information. 

Suppose, now, that we window xi(t) into two, 20 sections, xia(t) and Xib(t). 
These two new signals will now consist of only one frequency each, f 75 oo and fi 25 o. 
respectively. If we now take the Fourier transform of xia(t) and xib(t), separately, we 
obtain X|a(/) and X|b(/) as shown in Figure 3.2. Results show that for 0 < t < 20 ps, only 
f 7 !ioo is present, while for 20 ps < t < 40 ps, only fi 25 o is present. Hence, we now have 
reintroduced some temporal information into our analysis by windowing the data. This, 
time localization via windowing, is the basic idea behind the short-time Fourier transform 
(STFT) which provides a surface as the time-frequency structure which represents a given 
signal in the time-frequency plane [13]. 
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Figure 3.1: Time traces (a) of X](t) showing the change from f 75 oo to fiiso at t = 20 
\is and (b) of X 2 (t) showing the addition of two single sinusoid signals X 75 oo(t) and 
Xi 25 a(t) over the same time span 18 to 22 ps. The power spectrums of Xi(t) and 
X 2 (t) are shown in (c ) and ( d ), respectively, as being the same even though the 
two signals are different All frequencies are given in Hertz. 



Figure 3.2: Power spectrums (a) of Xi.(t) which is the portion of X|(t) from 0 to 20 ps, 
and (b) of X|b(t) which is the portion of X](t) from 20 ps to 40 ps showing how 
temporal information can be reintroduced into the signal analysis. Frequencies are 
in Hertz. 
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Assume that a signal, x(t), is stationary in a window, g(t), centered at time x. The 
Fourier transform of the windowed signal x(t) g*(t-x) leads to the STFT expansion given 
by: 

STFT(t,/) = J x(t)g*(t-'C) e dt , (3.12) 

where * denotes the complex conjugate. The STFT is a function of two variables, x and/, 
which represent time and frequency, respectively. Therefore, the STFT function maps a 
function, x(t), of only one variable t, representing time, into a two-dimensional time- 
frequency plane. 

There is a dual interpretation of the STFT as shown in Figure 3.3 [7]. The first 
interpretation is that of “windowing the signal” as illustrated by the sliding window, g(t), 
in Figure 3.3. In this interpretation, a window size is selected around some center time, x. 


SIdng 

Winoo<wg(t) 



Figure 3 J: The tiine>frequeiicy plane of the Short-Time Fourier Transform. It can be viewed 
either as a succession of Fourier Transforms of windowed segments of the signal (vertical bars) or 
as a modulated filter bank (horizontal bars). From Ref. [7] 
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and all frequencies within this window are computed. The second interpretation is that of 
a modulated window function, or a “filter bank” as illustrated by the implied modulated 
filter bank of Figure 3.3. Under this interpretation a bandpass filter is centered on some 
frequency, f, and all times within this filter are computed. With this dual interpretation in 
mind, let us call the bandwidth of the bandpass filter Af and the window, g(t), size At. 
Now, from the uncertainty principle we recall that the time-bandwidth product. At Af, has 
a lower bound given by; 

ArA/>^. (3.13) 

Ak 

Therefore, one cannot simultaneously obtain good time and frequency resolution, but has 
to trade one for the other. One potential drawback of the STFT is that this time frequency 
partitioning is flxed for the entire transform, as shown in Figure 3.4. Figure 3.4b shows 



Figure 3.4: Tiling of the time-frequency plane by the Short-Time Fourier Transform ( a ) a 
generic case and ( b ) the pattern from the example in Figure 2.2. 
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how the tiles are filled in the case of the example given in section E-B. E our signal of 
interest was composed entirely of transients or entirely of sinusoidal terms we would be 
able to choose a short window in the former case or a narrow filter in the latter case, and 
the STFT would prove to be adequate for our analysis. In practice, however, many 
signals have both long duration sinusoidal and transient components. For these types of 
signals the STFT is often not the right tool for adequate analysis. What is needed for 
these types of signals is a multiresolution analysis, one type of which may be obtained 
with wavelets. 


B. WAVELET ANALYSIS 

The Fourier series analysis equation, given in Equation 3.4, shows that the 
coefficients are formed by projecting the signal x(t) onto a set of complex exponential 

basis functions, e , where k = 0, ±1, ±2. One might call this form of analysis, 

“wave” analysis, since complex exponentials are wave functions (i.e., sines and cosines). 
In “wavelet” analysis, instead of using wave functions as our basis, we use wavelets. A 
wavelet, literally meaning “small wave” [14], has its energy concentrated in time [11]. 
Therefore, although it still has oscillating, wave-like, characteristics facilitating frequency 
analysis, its localization in time allows transient analysis [II]. 
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1. Wavelet Series 


This section is titled “wavelet series,” even though by convention the techniques 
about to be described are called wavelet transforms, to follow the analogy to Fourier 
analysis. We will remain consistent with the Fourier analysis terminology, but will also 
continue to point out conventional names found in the wavelet literature. This method 
will hopefully allow the reader to more closely follow the analogies, while still allowing 
comparisons to other discussions presented in the literature. The reader is encouraged to 
refer to Table 3.1 to avoid confusion. 

a. Continuous-Time Representation 

This representation of the wavelet series, called the discrete wavelet 
transform, will be developed by analogy with its counterpart in Fourier analysis, the 
Fourier Series. Consider the Fourier series equations pair reproduced here for 
convenience: 


Common 

name 

Consistent 

name 

Time, 

C or D 

Transform, 

C or D 

Input 

periodic 

Output 

periodic 

FS 

CTFS 

C 

D 

Yes 

No 

DFT 

DTFS 

D 

D 

Yes 

Yes 

DTFT 

DTFT 

D 

C 

No 

Yes 

FT 

CTFT 

C 

C 

No 

No 

DWT 

CTWS 

C 

D 

Yes or No 

Yes or No 

DTWT 

DTWS 

D 

D 

Yes or No 

Yes or No 

- 

DTWT 

D 

C 

No 

No 

CWT 

CTWT 

C 

C 

No 

No 


Table 3.1: ( C ), Continuous, and ( D ), discrete, periodic and non-periodic 
input and output relations for the Fourier and wavelet transforms. After 
Ref. [11] 




x(0 = (3.14) 

k=-oo 

a, = — dt (3.15) 

^ TO 


In the case of the Fourier series analysis equation, we have a set of basis functions, 
^-jkiKfoi ^ where ^ = 0, ±1, ±2, .... In wavelet analysis, we have a set of basis functions 
derived from our wavelet function, V|/(t), which is also referred to as the generating 
wavelet, mother wavelet, or analyzing wavelet [11]. As discussed at the beginning of 
section III-B, wavelets are useful for transient, as well as, dynamic sinusoidal signal 
analysis. This suggests that a wavelet transform will be a function of two parameters, 
namely frequency and time, much like the STFT. The first parameter introduced will 
position the wavelet in time by integer translations as shown below: 

WAt) = V/(t-r), t = 0,±1,±2. (3.16) 

The remaining parameter to be considered is frequency. Multiresolution analysis 
techniques allow the window size, At, and the bandwidth, Af, to vary across the time- 
frequency plane within the confines of the uncertainty principle. Variations are obtained 
in wavelet analysis by forcing the ratio of the bandwidth, Af, to center frequency, fc. to be 
constant. This constraint results in “constant-Q” filtering and the time-frequency plane 
partitioning shown in Figure 3.5 [7]. Figure 3.5 also illustrates the good time resolution, 
but poor frequency resolution obtained at high frequencies. Conversely, low frequencies 
display good frequency, resolution, but poor time resolution. 
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Figure 3.5: The tiling of the time-frequency plane by the wavelet transform. 

Thus, the time-frequency tiling described above is implemented by scaling the mother 
wavelet by a scaling factor, a, where a is inversely proportional to frequency, /. Our 
wavelet function can now be written as: 

where the xj-Ja term is used for energy normalization [7]. The relationships between 
scale, a, frequency,/, time resolution and frequency resolution are summarized in Table 
3.2. Now. restricting a to dyadic scales only (i.e., a-2') and let x = n2\ we can write our 


SCALE 

FREQUENCY 

FREQUENCY 

RESOLUTION 

TIME RESOLUTION 

low 


poor 

good 

mssil 

low 


poor 


Table 3.2: Relationships between scale, frequency and resolution. 
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wavelet series expansion of a given signal, x(t), as: 






(3.18) 


The two-dimensional set of coefficients, c,>, are called the discrete wavelet transform of 
x(t), and are analogous to the Fourier series coefficients, fl„, in Equation 3.4. 

We mentioned earlier that the wavelet function is also sometimes referred 
to as the ’’mother” wavelet. This term hints to the existence of another type of wavelet, 
the “father” wavelet. This function is more commonly referred to as the scaling function. 
We define the scaling function, (p(t), much as we did the wavelet function in Equation 
3.24, namely: 


or, again using a = 2‘ and x = w2' we have: 

Define the subspace of spanned by these scaling functions as: 

V,.Spa„j^.(i,]}. 


(3.19) 


(3.20) 


(3.21) 


Similarly, the subspace of spanned by the wavelet functions can be defined as: 



where Wi is the orthogonal complement of V/ in subspace V/.y. Therefore, the 
relationships of these subspaces is given by: 

=V^ © ®- • -©^2 © W;, (3.23) 

where is the final space spanned by the scaling function which is usually chosen to 
represent the coarsest detail of interest. The relationships between these subspaces are 
illustrated in Figures 3.6 and 3.8. Equation 3.23 allows us to rewrite Equation 3.18 in 
terms of both the wavelet and scaling functions for any function x(t) e l} as: 

+00 N 

X{t) = 2 + 2 (3.24) 

n = -oo i=\ /!=;-«> 

Note here that a(n) are the coefficients of the scaling function, (p„(t), which spans V„ of 
Equation 3.23. Similarly, d(i,n) represent the N sets of coefficients of the wavelet 
functions. which span subspaces, W), where i=l,2,...,N. 

The coefficients a(n) and d(i,n) define the discrete wavelet transform of 
x(t) in terms of both the scaling and wavelet functions, just as the coefficients, c,>, of 

V :d V Z> V =)V 3 



Figure 3.6: Orthogonal Vector Spaces V. and representing the spaces spanned by the scaling 
function and respectively. After Ref. [11] 
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Equation 3.18 did in terms of just the wavelet function. If these basis functions (p(t) and 
\ff(t) form an orthonormal basis or at least a tight frame [11], the coefficients may be 
derived as: 


a{n) = \x{t)(p„{t)dt, 

and: 


(3.25) 




(3.26) 


b. Discrete-Time Representation 

The nesting of subspaces, as shown in Figure 3.6, is achieved with the 
dilation equation given by: 

(p{t) = ^h(n)yf2(p{2t-n), /i = 0,±1, ±2,... (3.27) 

n 

where h(n) is the scaling function coefficients and maintains the norm of the scaling 
function as each successive dyadic scale is calculated [11]. Similarly, we can define the 
wavelet function in terms of the scaling function by: 

y/{t) = ^h^(n)^I2(p{2t-n), « = 0J:1, ±2,... (3.28) 

n 

where, 

A,(n) = (-l)"M/V-l-/i), (3.29) 

where N is the finite even length of h(n), as shown in Burrus [11]. Equation 3.27 can be 
interpreted as obtaining lower resolution by down-sampling by two after lowpass filtering 
with the half-band filter with impulse response h(t). Similarly, Equation 3.28 can be 
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interpreted as obtaining lower resolution by down-sampling by two after highpass 
filtering with the half-band filter with impulse response hi(t). The analysis filter bank just 
described are illustrated in Figure 3.7. 

The coefficients of Equations 3.27 and 3.28, namely h(n) and hi(n), can be 
viewed as digital filters and the coefficients of Equation 3.24, namely a(n) and di(n), can 
be viewed as digital signals. With this view it can be shown [11] that the coefficients of 
Equations 3.25 and 3.26 become: 

«,+i (n) = X ” 2n)fl, (m), (3.30) 

m 

and: 

di^.^{n)=^^h^(m-2n)al(m), (3.31) 

m 

respectively. The above discussion leads to what is commonly referred to as the filter 
bank view of the DWT or the discrete time wavelet transform (DTWT) [11]. These last 
two equations show that we obtain the coefficients at the next level of scale, i+I, by 
convolving the approximation coefficients, ai(m), at level i with the time reversed filter 
coefficients. h(-n) and hi(-n), and then down-sampling by two. This idea is illustrated in 
block diagram form in Figure 3.7, where -12 indicates decimation, or down-sampling, by 
two and the other blocks represent filtering operations. As mentioned previously, h(n) is a 
lowpass filter and hfn) is a highpass filter. Originally, we have a signal, xfnj. For this 
case, if the samples of the input function, xfnj, are above the Nyquist rate, they represent 
a good approximation to the scaling coefficients at that scale, meaning that the 


28 






coefficients, oq, equal x[n] [11]. Next, we filter ao with h(n) and hj(n) and down-sample 
by two, to get the next scale level coefficients, and di, respectively. The resulting 
coefficients represent the next coarser scale level which give better frequency resolution, 
but poorer time resolution. We can repeat the process with the approximation 
coefficients, ai, to get 02 and d 2 , as indicated in Figure 3.7, which represent the next 
coarser scale level with even better frequency resolution at the price of even poorer time 
resolution. Normally, this process is repeated until the desired level of coarseness, or 
frequency resolution, is obtained. 

Wavelet transforms can also be described in terms of filter banks as shown 
in Figure 3.7. Let i=0 for the example to follow. When the sampled signal is passed 
through the first filter bank consisting of hi(-k) and h(-k), the highpass filter and lowpass 
filter, respectively, the original space is divided into Wj and Vj as shown in Figure 3.8a. 
After passing through the second stage of the filter bank, Vj is divided into W 2 and V 2 . If 
we assume that we have now reached our coarsest scale desired, we now have the entire 
signal space of interest represented by W, 0 Wj ® Vj, as shown in Figure 3.8b. If we have 



Figure 3.7: Analysis filter bank for calculating the discrete time wavelet transform. 
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2 

(a) 


\H(o))\ 
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(b) 



8 4 2 

(c) 

Figure 3.8: Spectral Partitioning by analysis Hlter bank. 


16 8 4 2 

(d) 


not reached the coarsest scale desired, we continue to iterate the filter bank to further 
partition the frequency bands, as shown in Figures 3.8c and 3.8d. Note that the gain of 
the frequency response in Figure 3.8 has been normalized for all bands. In actuality, the 
areas under the curves would be equal in order to meet the requirements for conservation 
of energy. 


Synthesis is performed by simply reversing the process as shown in the 
filter bank representation of Figure 3.9 and given by: 

a, in) = “ 2m) + “ 2m). (3.32) 

IK m 
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Note in Figure 3.9 that T 2 indicates up-sampling by two by inserting a zero between 
each point, and that the result is filtered with h(n) and h](n) not their time reversed 
counterparts. The analysis Equations 3.30 and 3.31 taken together with the synthesis 
Equation 3.32 are known as Mallat’s algorithm [15]. Recall that the two transforms we 
have just introduced, called the DWT and the DTWT, are in the Fourier analysis 
terminology called the continuous time wavelet series and the discrete time wavelet 
series, respectively. Therefore, we are left to draw an analogy with the “continuous” 
Fourier transforms (see Table 3.1). 

2. Wavelet Transform 

tL Continuous-Time Representation 

Recall that the STFT of Equation 3.12 is of the form; 

STFT{r,f)= rxiOyf * it-r)dt, (3.33) 

where 1^(0 = g(t)e~^'’^'^. If we introduce the scaling factor a into Equation 3.33, we 
obtain the continuous time wavelet transform (CTWT) expansion given by: 



Figure 3.9: Synthesis filter bank for calculating the discrete time wavelet transform. 
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1 f+“ ft — T 

CWT(T,a) = - 7 =J xit)iff*\ - 


\dt, 


(3.34) 


where a = f^j f and the term l/ 4a term is introduced for energy normalization [7]. 

Although not used in this work, one can define the inverse CWT as: 


x{t) = K\\\cWY{a,T)\if\ 
a 


V a 


dadx , 


(3.35) 


where the normalizing constant is: 



and T(/) is the Fourier transform of y/(t) [11], 


(3.36) 


b. Discrete-Time Representation 

Signal processing in most real world applications is performed on a digital 
computer. This implies that we must use the discrete, sampled, version, x[nj, of the 
continuous signal, x(t), which we wish to process. This necessity leads us to the need for 
a discrete time continuous wavelet transform (DTCWT), or to be consistent with Fourier 
analysis terminology, simply a discrete time wavelet transform (see Table 3.1). The 
actual analysis and synthesis Equations for the DTCWT will not be presented here, but 
rather the important conceptual differences with the DWT will be discussed. 

Recall that for the DWT that the scales, a, are restricted to dyadic scales 
(i.e., fl = 2', for / = 0, ±1, ±2,...) to allow for the filter bank implementation. However, 
no such restriction needs to be made for the DTCTWT. The scale ranges, continuously. 
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from that of the original signal, up to whatever maximum the user desires [16]. This 
continuous range of scale is the first important difference compared to the DWT. The 
second significant difference is that shifting of the wavelet function by the DTCWT is 
continuous [16]. The analyzing wavelet is shifted just as smoothly over the time domain 
of the signal being analyzed as it is scaled over the frequency domain. In the next chapter 
we will look at the abilities and limitations of the DTCWT in detecting transients by 
using the Wavelet Toolbox’s implementation of this transform. 

So far we have introduced frequency hopping signals and the various 
forms of wavelet analysis. In the next chapter we will apply the latter to the former to see 
how wavelet analysis might be used in the detection and estimation of frequency hopping 
signals. 
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IV. WAVELET BASED TRANSIENT DETECTION 


The frequency hopping detection scheme considered in this work involves the 
detection of discontinuities. This chapter, first, briefly reviews how wavelets can be used 
to detect signal discontinuities. Next, it defines the temporal correlation function used in 
our detection scheme. Finally, it introduces several preprocessing tools used to minimize 
potential performance degradation when additive white Gaussian noise is present in the 
communication signal. 

A. DETECTING DISCONTINUITIES USING WAVELETS 

Consider the three signals given in Figure 4.1. Figure 4.1a represents a signal, 
x(t), with a discontinuity. Figure 4.1b represents a signal with a discontinuity in the first 
derivative, and Figure 4.1c represents a signal with a discontinuity in the second 
derivative. Wavelets may be used to detect each of these discontinuities, if the chosen 
wavelet is able to represent the highest order derivative present in the signal function, as 
any wavelet with, at least, p vanishing moments can be used to detect a discontinuity in 
the p-\ derivative [10]. However, the wavelet’s ability to detect these discontinuities 
erodes quickly in the presence of noise, as will be shown later. 
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Figure 4.1: Function discontinuities (a ) in the function (step function), ( b) in the first derivative 
(interrupted ramp function), and ( c ) in the second derivative (frequency hopping function). 

Now, consider the various wavelets shown in Figure 4.2. The Morlet wavelet is an 
example of a continuous wavelet that can represent all three of the functions in Figure 
4.1, if appropriately compressed or dilated (i.e., it can be compressed to until it contains a 
vertical portion, dilated to match a ramp region, or dilated less to match the sinusoidal 
terms). The Haar wavelet shown in the figure has only one, p=l, vanishing moment and, 
therefore, can only detect discontinuities in the zero, p-l=0, derivative (i.e. the function 
itselO- Further, note the Haar wavelet does contain a portion which can match functions 
with discontinuities, such as the step function of Figure 4.1, exactly, but no portions 
which can match the ramp or sinusoidal functions of Figure 4.1. Therefore, it is not able 
to detect discontinuities in either the first or second derivatives, as shown in the next 
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Moliet 


Haar 



Figure 4,2: The ( a) Morlet, ( b ) Haar, ( c ) Daubechies-4 and (d ) Daubechies-6 wavelets with 
their scaling functions shown as dashed lines. 

section. Similar arguments hold for the Daubechies-4 (p=2) and Daubechies-6 (p=3) 
wavelets shown in Figure 4.2 which are able to detect discontinuities in functions or their 
derivatives up to the first and second degree, respectively. 

1. The Continuous Wavelet Transform 

As mentioned earlier, the continuous wavelet transform (CWT) can be used to 
detect discontinuities in a function and its first and second derivatives. The discrete time 
continuous wavelet transform of Subsection III.B.2.b was used to detect the 
discontinuities of the functions shown in Figure 4.1, using the Morlet wavelet. Figure 
4.3. shows that it is able to detect discontinuities in the function shown in Figure 4.1a, in 
the first derivative of the function shown in Figure 4.1b, and in the second derivative of 





















the function shown in 4.1c. However, since the CWT is computationally expensive, we 
will choose not to use it in our detection algorithm. 

2. The Discrete Wavelet Transform 
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Figure 4J: The first level detail coefficients, di, from the Continuous Wavelet Transform of the 
functions of Figure 4.1, ( a ) the step function, ( b ) the interrupted ramp function, and ( c ) the 
frequency hopping function, using the Morlet wavelet of Figure 4.2a. 



Recall that the Daubechies-6 wavelet has three vanishing moments allowing it to 
detect discontinuities in a function and discontinuities in both the first and second 
derivative of a function. Figure 4.4 shows that the Daubechies-6 was indeed able to 
detect such discontinuities in the functions shown in Figure 4.1 using the DTWT of 
section m.B.I.e. Further, recall that the Daubechies-4 and the Haar wavelets have only 
two and one vanishing moments, respectively. As shown in Figure 4.5, the Daubechies-4 
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wavelet was only able to detect the discontinuities in the step and interrupted ramp 
functions, but not the discontinuity in the second derivative of the frequency hopping 
function. Figure 4.6 shows the results of performing a DWT on the functions shown in 
Figure 4.1 using a Haar wavelet. As expected, it was only able to detect the discontinuity 
in the step function. 
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(C) 

Figure 4.4: The first level detail coefTicients, d|, from the Discrete Wavelet Transform of the 
functions from Figure 4.1, ( a ) the step function, ( b ) the interrupted ramp function, and ( c ) the 
frequency hopping function, using the Daubechies-6 wavelet of Figure 4.2d. 
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Figure 4.5: The first level detail coefficients, di, from the Discrete Wavelet Transform of the 
functions of Figure 4.1, ( a ) the step function, ( b ) the interrupted ramp function, and ( c ) the 
frequency hopping function, using the Daubechies-4 wavelet of Figure 4.2c. 
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Figure 4.6: The first level detail coefficients from the Discrete Wavelet Transform of the 
functions in Figure 4.1, ( a ) the step function, ( b ) the interrupted ramp function, and ( c ) the 
frequency hopping function using the Haar wavelet of Figure 4.2b. 
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3. The Continuous Wavelet Transform Applied to Functions with Additive 
Noise 

Adding white Gaussian noise (AWGN) to the signal, x(t), to obtain a SNR of 60 
dB, results in the signal shown in Figure 4.7. Note that the noise is not even perceptible 
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Figure 4.7: Functions of Figure 4.1 embedded in 60 dB of AWGN. 


with the given figure’s resolution. However, by examining the wavelet transforms of 
these functions with this moderate noise, insight may be gained into the robustness of 
these discontinuity detection techniques. Figure 4.8 shows the results obtained when 
using the DTCWT with the Morlet wavelet. Note, that even though it is still able to 
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Figure 4.8: The first level detail coefficients, dj, from the DTCWT on the functions of 
Figure 4.7 using the Morlet wavelet. 


detect a discontinuity in the step function in noise, it fails to detect discontinuities caused 
by the first and second derivatives of the functions of Figures 4.7b and 4.7c, respectively. 

4. The Discrete Wavelet Transform Applied to Functions with AWGN 

Applying the DTWT using Daubechies-6, Daubechies-4, and Haar wavelets to 
the functions shown in Figure 4.7 leads to the results shown in Figure 4.9 to 4.11. We 
note that discontinuities in the first and .second derivatives are not detected due to noise 
degradations present in the signals. Therefore, one can conclude that a scheme based on 
detecting discontinuities in the function itself, vice its derivatives, is more robust. 

Thus, the next question becomes: when all these wavelets can detect these 
di.scontinuities in the function itself, how does one decide which wavelet to use when 
looking for such discontinuities? The DTWT computed using filter coefficients, h(n) and 
hi(n), of short length are not only less expensive computationally, but also produce better 
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time localization [10]. Based on the latter consideration the Haar wavelet should be a 
good candidate, since it has the shortest filters of the three. 
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Figure 4.9: The first level detail coefficients, d|, from the DTWT of the functions of Figure 4.7 
using the Daubechies-6 wavelet 


5. Averaging of Scales 

The averaging of several scales can enhance a wavelet’s ability to detect 
discontinuities in noise. The idea is that for true discontinuities, the spikes will line up 
across all scales, while the spikes due to noise will not line up. Consider, for example, the 
stair step function shown in Figure 4.12a and the same function in noise in Figure 4.12b. 
The Haar wavelet detects the discontinuities perfectly in the no-noise case, as shown in 
Figure 4.12c. However, the performance is seriously degraded by the noi.se. as shown in 
Figure 4.12d. 
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(c) 

Figure 4.10: The first level detail coefticients, di, from the DTWT of the functions 
of Figure 4.7 using the Daubechles-4 wavelet. 
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Figure 4.11: The first level detail coefficients, d|, of the DTWT of the functions of Figure 4.7 
using the Haar wavelet. 
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Figure 4.13 shows the DTWT of the stair step function in noise at scales 1 through 
6 computed using the Haar wavelet. Note that all step times can be detected if the 
wavelet coefficients obtained from the first six scales are summed and clipped at some 
threshold, as shown in the top middle plot of Figure 4.13 labeled “Output of detector” 
[ 10 ]. 


Input Signal 


Input Signal 



Haar Wavelet, level 1 


Haar Wavelet, level 1 


(a) 


(b) 


(c ) 


Figure 4.12: Detection of the stair step function (a ) with no noise, ( b ) in noise with the 
corresponding detector output shown in (c) and (d ), respectively. From Ref. (10]. 
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Figure 4.13: Detection of stair step function in noise using scale averaging. From 
Ref. [10]. 
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B. 


THE TEMPORAL CORRELATION FUNCTION 


The temporal correlation function (TCP) of a signal x(t) is defined as: 

f f 

TCFrt,T) = A: t + - ;c* t-- , (4.1) 

y 2j y 2) 

where t is the absolute center time and t is the lag time. This function has been shown 
useful when analyzing nonstationary signals [17]. Note that the TCP is conjugate 
symmetric along the T axis as: 

TCp(t, t) = TCP* (t -t) , (4.2) 

where * represents the complex conjugate [17]. Therefore, no additional information is 
gained by calculating the function for negative values of the lag, T. In addition, note that 
the TCP of a real signal may exhibit interference terms due to its nonlinear definition. 
Por example, consider a real sinusoidal signal x(l): 

x{t) = sin(2;r f t +6). (4.3) 

The lesulting TCP of x(t) is given by: 

TCF(t.T) = ^[cos(2;r/t)-cos(2;r[2/]t + 20)]. (4.4) 

The first term inside the square brackets consists of auto-terms at the frequency,/, of x(i), 
while the second term consists of interference terms at twice the frequency,/, of x(t) [17]. 
Note, however, that the TCP of the analytic signal, xjt), obtained from x(t} given by: 
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= (4.5) 

leads to the TCF function: 

TCF, (t,T) = e'''^^ (4.6) 

Equation 4.6 shows that the TCF expression contains only auto terms at the frequency of 
xjt) [Fried]. For this reason the analytic form of the signal, x(t), is often preferable to its 
real counterpart and is the form we will choose to work with in our detection scheme to 
be developed. 

Next, consider the nonstationary analytic signal given by: 

x,(t) = e‘^’'"[u(t)-uU-T^)]+el‘l''[u(t-T^+\)-u{t-T)], (4.7) 


for 0 < r < r, where Thop is the time of the hop (or change in frequency) from// to/ 2 , and 
where u(t) is the unit step function given by: 


f 1 ,for r > 0 

[O , for t < 0 


(4.8) 


Substituting Equation 4.7 into Equation 4.1 leads to: 
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7CF(t,T) = 
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+ i'^A 
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(f2-hY*[^~y 


r+--7 +1 

V 2 


V 2. 


V 2 


hap 


+ul /+‘~-r 

2 y 


V 2 *‘’V 


rCF,(t,T) + 7’CF2(t,T)+ 7CF,2(t,T), 
where TCFi(t, T), TCF 2 (t, t), and TCFi 2 (t, t) represent the 1®', 2"**, and terms of 
Equation 4.9. Figure 4.14a presents a phase plot of Equation 4.9 for some arbitrary // and 
f:. The combinations of the different shifted versions of the unit step functions, u(t), force 
the TCF to only take on values within the regions shown in Figures 4.14b to 4.14d. Note 
that TCF|(t, T) is a function of// and Tonly. The second term, TCF 2 (t, T), is a function of 
/2 and Tonly, while the last term, TCFi 2 (t, T), is a function of//,/ 2 , t, and t. Note that the 
frequency hopping time, Thop = 100, is located where the region covered by TCFi(t, T) 
ends and the region covered by TCF 2 (t, t) begins. Further, note that for a given value of 
T. the terms within the triangles (i.e., the regions where TCFi(t, T) and TCF 2 (t, T) are 
defined) are con.stant, although at different levels. This fact is further illustrated in Figure 
4.15a for a few values of t, and will be exploited in our detection algorithm. Equation 
4.9 showed that the phase behavior within the cross-terms region, TCFi 2 (t, T) is linear. 
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The phase expression is periodic over 2n producing discontinuities at regular intervals 
across TCFi 2 (t, t). It is important to realize, however, that the period of these intervals is 
a function of// and /2 and, therefore, not predictable without knowing/; and/ 2 , which in 
general we do not. Nevertheless, for any given value of the lag, T < Thop, this region of 
cross-terms is centered on the hop time, Thop, another fact which may be exploited. 

It should be further highlighted that within a given auto-term triangle, there is only 




(c) (d) 


Figure 4.14: Phase plot of the Temporal Correlation Function obtained from x,(t) defined in 
Equation 4.7. ( b ) shows the portion of ( a ) which contains the auto*terms, TCF|, due to// only, ( c ) 
shows the portion of (a ) which contains the auto-terms, TCF 2 , due to/; only, and ( d ) shows the 
portion of ( a ) which contains the cross-terms, TCFu, due to both// and//. 
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one frequency component along the t axis. Therefore, if one were to calculate the 
frequency response within one of these triangles, across the values of t, an estimate of 
the frequency during the period of t, covered by the given triangle could be extracted. 

So far, only signals without the presence of noise have been considered. In 
practice, some amount of noise is always present and the following section introduces 
additional preprocessing tools designed to facilitate the extraction of frequency hopping 
features in noisy environments. 

C. PREPROCESSING TECHNIQUES 

As mentioned before, noise added to the signal impedes our ability to detect 

discontinuities in frequency hopping signals, or in general, to detect transients. This 
section investigates how the additive noise alters the TCP expression, and presents 
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Figure 4.15: Phase behavior of the TCF at lags 16 (top), 48 (middle), and 60 (bottom) for ( a ) the 
case with no noise and ( b ) the case with additive WGN. 
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several preprocessing techniques to minimize the noise effects. Specifically, we will 
investigate how the noise affects the phase information in the TCF computed on the 
analytic signal, Xa(t), given in Equation 4.7. Recall that the frequency value information 
and the hopping time may be obtained using the phase of the TCF computed from an 
analytical frequency hopping signal, x(t}, as shown earlier. In addition, note that: 

1) the phase components of the auto-terms, TCFi(t,T) and TCF 2 (t,T), have 
constant values for a given lag time v, 

2) the phase components of the cross-terms, TCFi 2 (t,T), are linear in t for a given 
lag time t. 

Figure 4.15 plots the phase information obtained from the TCF expression shown in 
Figure 4.14 for lags T= 16, 48 and 60 in noise free and noisy environments. Note that the 
hopping time {Th„p = 100) occurs in the middle of the linear phase excursion. Also, note 
that the noi.se causes random spikes in the phase information. Such spike occurrences 
must be minimized if we are to produce a robust, reliable detection algorithm. 

1. Unwrapping the Phase 

The pha.se, p(t), of a signal, x(t), may be unwrapped as: 

P(0 if |p(/)-p(r-l)|<;r, 

unwrap(p(/)) = jp(/) + 2;r if p{t)-p{t-\)<-n, (410) 

p{t)-2n if p(t)-p(i-\)>7:. 
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Unwrapping the phase serves two purposes. First, it removes the discontinuities present 
in the region covered by TCFi 2 (t,T). The resulting unwrapped phase becomes purely 
linear for a given lag time t in the region covered by TCFi 2 (t,T), as shown in Figure 
4.16a. Second, the unwrap function is, also, very effective at removing the random spikes 
in the phase due to noise, as shown in Figure 4.16b. However, as we shall see, the 
wavelet transform is sensitive to even low level noisy behavior. The same properties 
which make the WT adept at detecting transients make it susceptible to false alarms due 
to the noise. Therefore, further “noise quieting” techniques prove beneficial. 

2. Median Filter 

A median filter of size N, is given by: 
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Figure 4.16: Unwrapped phase behavior at lags requal to 16 (top), 48 (middle), and 60 (bottom) 
for (a) the case with no noise and (b) the case with additive WGN (SNR = 10 dB). 
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(t) = median(x(OwAr (0) 


(4.11) 


where w//t) is a window of length N centered on t. The median filter is a nonlinear filter 
with the important characteristic that it does not average and, therefore, preserves 
discontinuities. Examples of the effects of the median filter on our noisy phase 
information will be given in the next chapter. 


3. Differentiation 

Consider a function which is at a constant value, Ci, for some period of time, 
0<t<t,, then ramps up with slope, m, for another period, t, < t , and is then 
assumes another constant value, C 2 , for the remainder of its duration, t 2 <t <T, as 
shown in Figure 4.17a for Ci= -1.2474, € 2 = 13.832, t/=70, ^2=130, and 7’=256. For future 
reference, note that this figure closely resembles that of Figure 4.16a (bottom) with the 
exception that the end terms, which will be discussed later, have been set to the same 
constant value as those before and after the ramp (i.e. the constant levels due to the auto¬ 
terms of Figure 4.14a). Differentiating the function shown in Figure 4.17a, leads to 
Figure 4.17b as the constant terms go to zero, and the ramp portion becomes a constant 


C, - C 

equal to the slope, m = — -^ = 0.251, of the ramp. Therefore we have created a pulse. 


C -t. 


a function with two discontinuities, of height m = 0.251 and width equal to t:-1/ = 60 as 
shown in Figure 4.6b. This is another key feature which we will exploit in our detection 
algorithm. 
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The corresponding desired results obtained from the application of these functions 
are listed as follows: 

1. The TCF produces a two-dimensional correlation function with the property 
that when applied to a analytic signal such as the one in Equation 4.7 (i.e. a 
complex sinusoid with one frequency hop) it results in a pattern, as shown in 
Figure 4.14. Furthermore, the cross-terms, TCF 12 , are centered on the hopping 
time, i=Thop, for 0< t < T and 0<r<T^p. 

2. When the TCF phase is unwrapped, for any given lag, t, it results in a signal 
which is constant for the period of the first auto-terms, TCFi, then is a ramp 
for the period of the cross-terms, TCF 12 , and then is constant again, at a 
different level, for the period of the second auto-terms, TCF 2 . 



(a) 


S 0.2 
« 0.1 
£ 0 
- 0.1 


100 200 
Time (samples) 


( b ) 


Figure 4.17: A ramp function is shown in ( a ), and its 
derivative in ( b ). 
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3. The application of a median filter to the phase information in noise can 
smooth out the noise, but with the important characteristic that it will preserve 
the discontinuities. 

4. Differentiation of a function of the type described in item two above results 
in a pulse centered on the hopping time, Thop, for t< Thop- 

Having explored the abilities and limitations of using wavelet transforms to detect 
discontinuities, and having examined some preprocessing tools, we are now ready to 
develop our frequency hopping signal detection algorithm. 
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V. DETECTION AND ESTIMATION ALGORITHM AND SIMULATION 

A. DETECTION AND ESTIMATION ALGORITHM 

Using the tools described in the last chapter, the algorithm steps for the detection 
and estimation of frequency hopping signals in noise can now be enumerated. The steps 
are as follows: 

1. Transform the real signal into an analytic signal. 

2. Segment data into frames of length less than or equal to the minimum 
hopping time, Thop_min- This assumption ensures that, at most, one hop will 
be present in the processing frame. 

3. Compute the temporal correlation function on each frame. 

4. Extract the phase information by calculating the angle of the TCP. 

5. Unwrap the phase of the TCP along the time axis, t. This step is done to 
remove artificial discontinuities in the cross-terms region due to the phase 
being periodic in 2jt. Unwrapping the phase, also, reduces the apparent 
random spikes in the phase due to noise. 

6. Apply a median filter to the phase of the TCP along the time axis, t, of length 
five. This step is done to reduce the noise effects prior to differentiating, 
since differentiating accentuates these effects due to noise. 

7. Differentiate the phase information along the time axis, t. This step changes 
the unwrapped phase of the TCP from a ramp function to a pseudo-pulse. 
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8. Apply another median filter along the time axis, t, of length 25. The length 
of 25 has proven to work well with the Thop_mn chosen for our simulations 
which will be described below. This step is done to again remove the effects 
due to noise which were accentuated by the differentiation operation in step 
seven. 

9. Calculate a discrete wavelet transform along the time axis, t (i.e. for each lag, 
T, of the TCP) using the Haar wavelet. This step is done to detect the 
discontinuities at the edges of the cross-terms region. 

10. Sum the wavelet coefficients of the first two scales, di and da of the DWT. 
Since each successive scale is down-sampled by two when performing a 
DWT, this step is accomplished by setting d2(2n+l) = d2(2n) for n = 
0,1,2,...,N-1 where N = lengthCda), before performing the summation. 

11. Perform a 457135° summation across all values of lag, x, to obtain a 
detection vector which has time as its index. This step will be explained 
further in Subsection V.C.6. 

12. Threshold the data in the resulting detection vector obtained in step 11. 
Further details regarding this threshold scheme are presented in Appendix B 
and Subsection V.C.7. 

13. If the thresholded detection vector contains peaks, the maximum peak value 
time index represents the estimated hopping time. 
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With the detection and estimation algorithm given above, we will use the remainder of 
this chapter to describe the simulation process, apply the algorithm to an example frame, 
and, finally, to provide the simulation results. 

B. SIMULATION 

Simulations were conducted to test the effectiveness of the detection and 
estimation algorithm given above. The Matlab source code for conducting the 
simulations is given in Appendix A. Five hundred trial experiments were conducted in 
six different signal to noise ratios (SNR) ranging from 15 dB down to -3 dB. The basic 
idea behind the experiments was to simulate signals that had already been segmented, as 
specified in steps one and two of the algorithm. The problem then becomes to determine: 
a) whether or not a frequency hop exists within the given frame; b) to estimate the 
hopping time when a frequency hop is detected. 

1. Assumptions 

As stated in the introduction, it was desired to make as few assumptions as 
possible on the nature of the frequency hopping signal. With this goal in mind the 
a.ssumptions were limited to three. The first assumption is of a frequency range within 
which the spread spiectrum signal remains. For the simulations conducted, this range was 
assumed to be from 1 MHz to 24 MHz. The second assumption was of a minimum 
hopping time, Thop.mm. which was chosen to be 256 sample points. At a sampling rate of 
50 MHz, this translates into a minimum hopping time of 5.12 ps. The third, and final. 
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assumption was on a minimum frequency differential, Af, for the hop which was chosen 
to be 1 KHz. 


2. Signal Generation 

The signals were generated by first choosing a random hopping time, Thop, 
between sample point 26 and 231 of the 256 point sample frame. The first and last 25 
points were disallowed as candidates due to problems with edge effects. In practice, an 
overlapping scheme could be used for full coverage of the signal. If it is determined that 
a hop will occur within a simulation frame (i.e. Thop ^ 0), then both first and second 
frequencies (i.e. frequencies before and after the hop) are randomly generated. The result 
is a signal with, at most, one hop which can be from any frequency,/;, to any frequency, 

fz, such that 1 MHz <f,,f 2 ^ 24 MHz, and |/, - /J > A/ . 

The SNR, in decibels, is defined in the simulations as: 


SNR = l0\og 


10 


^ Psignal ^ 
V ^ noise J 


(5.1) 


where o denotes the standard deviation squared, or variance, of the noise, and Psi^nai is the 
signal power. Forcing the signal to be of unit amplitude. Equation 5.1 becomes: 


SNR = lOlog 


10 


1/2 

y^naise J 


(5.2) 


which allows the calculation of the standard deviation of the additive white Gaussian 


noise as: 
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rj -SNR 

a . =J--10 >0 . (5.3) 

''noise V 2 ' ^ 

C. APPLICATION OF DETECTION AND ESTIMATION ALGORITHM TO 

A SAMPLE FRAME 

1. Choosing a Sample Frame 

We choose the eighth frame of the signal defined earlier in Section II-B, for the 
example. This frame is 256 sample points long from sample f793 to sample 2048, as the 
frequency hopping signal is segmented into frames of size Thop_min equal to 256 points. 
Therefore, the eighth frame starts at point 256-7 + 1 = 1793, and ends at point 
256 • 8 = 2048. The hopping time for this frame is at 40 |xs, which corresponds to point 
2000 for a sampling frequency of 50 MHz. Thus, the hopping time is located at time 
sample 2000 - 1792 = 208 inside the frame. The hopping time within our sample frame is 
called Tiuip. The frequencies on either side of the hop strefnso = 6.250 MHz and ft 7727 - 
22.727 MHz. Finally, white Gaussian noise is added to the signal to obtain a SNR equal 
to 10 dB. 

2. The Temporal Correlation Function 

Figure 5.1 plots the phase of the TCF expression for the analytical function 
derived from the communication signal. Recall that the bottom triangle is due to the 
auto-terms at frequency fi 25 o and the top triangle is due to the auto-terms at frequency 
f 17727 . while the rest are cross-terms due to both frequencies. Figure 5.1b shows the 
values of the TCF phase for a constant lag equal to 30, Tso- 
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3. Extracting Phase Information from the TCP 

Steps four and five of the detection and estimation algorithm manipulate the phase 
information of the TCP into a convenient form using the preprocessing tools introduced 
in Chapter fV. Figure 5.2 plots the unwrapped phase of the TCP shown earlier in Figure 
5.1. Figure 5.2b shows the effect of unwrapping the phase at lag, tso. 


Phai»o?TCFolxB[n] 



O SO 100 150 aOO 230 

T»n» funolpa) 


(b) 

Figure 5.1: The ( a ) phase of the Temporal Correlation Function computed on an analytic 
frequency hopping signal in noise and ( b ) a closer view at the constant value of lag equal to 30. 
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Figure 5.2: The unwrapped phase information of the TCF of Figure 5.1. 


4. Constructing the Pulse 

Steps six to eight in the detection and estimation algorithm transform the phase 
information into a pulse-like signal, or pseudo-pulse. Note that the phase information is 
usually very sensitive to noise degradation. Thus, several steps are added to minimize 
noise effects. We first apply a median filter of length 5 to minimize the effects due to 
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noise while preserving the global trend of the phase information. Next, we differentiate 
the phase along the time axis t to obtain the pseudo-pulse. Finally, we apply a second 
median filter of length 25 to again minimize noise effects while preserving the pulse-like 
shape of the function. Figure 5.3 plots the resulting signal obtained from the TCF phase 
for the lag time Tjo- We now have constructed a pseudo-pulse having the width of the 
cross-terms and, roughly, centered on Thop- A few comments can be made; 

• when the signal is noise free, such processing leads to a perfect pulse of a height 
equal to the slope of the unwrapped TCF phase, a width equal to the width of 
cross-terms, and centered on Thop (i-e. centered at sample point 208 of the time 
axis). 

• discontinuities at each edge of the pseudo-pulse are preserved when the signal is 
filtered with the median filter of length 25. 


|2 

0 

0 so 100 150 200 250 

Time (samples) 

Figure S3: The pseudo-pulse formed by differentiating and median filtering the signal of Figure 
5.2b. 







5. Detecting Discontinuities with the DWT 

Next, steps nine and ten, of our detection and estimation algorithm are applied to 
the pseudo-pulse function. Step nine computes the discrete wavelet transform using the 
Haar wavelet. The transform is well matched to detect discontinuities at the leading and 
trailing edges of the pseudo-pulse. Next, the wavelet coefficients obtained for the first 
two scales are averaged to minimize the noise degradation, and to enhance the probability 
of detecting the discontinuities, as shown earlier in Subsection V.D.4. Figure 5.4 plots 


DWT on Preprocassed TCP of x8In] 



DWT on P r apro c— od TCP at 30 



(b) 


Figure 5.4: DWT coeflicients computed for each value of T (a) on the preprocessed phase of the 
TCF(t,T), and ( b ) for x = 30. 
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the resulting wavelet transform obtained from the TCF phase information on the entire 
TCF, and a detailed view of the transform obtained for lag T 30 . 


6 . Constructing the Detection Vector 

Step 11 of our algorithm constructs what will be called the detection vector by 
performing a 45°/135° summation on the DWT matrix. The Matlab code for this 
summation can be found in Appendix A. However, the basic idea is to sum all the values 
which represent the edges of the cross-terms in the TCF at 45° and 135°, so that they 
reinforce each other only at Thop, which we recall is located at point 208 in this example. 
Figure 5.5a shows pictorially the effects of this summation. The arrows in Figure 5.5b 
show the directions of summation arrows pointing upward representing 45° and the 


DWT on R rop r oo —o ocf TC^ o# xSfnJ 



Tb« Dolootton V«ctor 



100 

Tim* (ft 


(C) 

Figure 5.5: Detection vector constructed by performing a 45®/135® summation across the columns 
of the signal shown in Figure 5.4a. ( a ) shows a representation of the effects of the summation 
shown in ( b ) where the arrows indicate the direction of summation. ( c ) shows a plot of the 
actual detection vector. 







arrows pointing downward representing 135°. The actual detection vector is plotted in 
Figure 5.5c. Note that the peak indeed appears to very close to point 208. 

7. Threshold Determination 

Once the detection vector has been formed a decision must be made as to whether 
or not a hop has occurred within the frame. The statistics of the empirical data from the 
experiments suggested that the variance of the detection vector would be the best 
indicator of whether or not a hop had occurred. As a result, the threshold, TthreshoU, is 
chosen as a multiple, k, of the variance of the detection vector when no hop has occurred 
within the frame as: 

= k ■ vai(delection vector„„_^^ (t)). (5.4) 

The threshold determination was also guided by the fact that the cost associated with the 
probability of a missed detection. Pm = [1 - probability of detection (Pp)] far outweighs 
the cost associated with the probability of a false alarm, Pfa, as the hopping time 
estimation is only the first step in a complete frequency hopping signal detection scheme. 
Note that once the hopping times are estimated, the signal frequencies need to be 
extracted to demodulate the actual message. This step can easily be done by applying 
frequency analysis to the estimated hopping intervals. Thus, in the case of false alarms, 
frequency analysis would show the same frequencies in two, or more, consecutive 
hopping intersals. resulting in no message degradation. However, a missed hopping time 
will result in degradations in the frequency estimation step, and errors in the decoded 
message. 
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Receiver operating characteristics (ROC) curves were generated for each of the 
six SNR levels and an appropriate threshold chosen for each. This curve, along with the 
plots of Pd vs. Tthreshoid and Pfa vs. T,hreshoid, is shown in Figure 5.6 for the SNR = 10 dB 
case. Similar plots for the other SNR levels are shown in Appendix B. When the SNR = 
10 dB, the threshold is chosen to be jt=30 times the variance of the detection vector 
generated in a “no hop” frame, which leads to “ 1-8907 x 10”*. Simulations show 

that the detection vector variance is equal to 0.0076, which exceeds the selected 
threshold, T,hreshoid = 1.8907 x 10'^. Therefore, a hop is detected. The value of the hop 
time is estimated next by locating the time index of the detection vector maximum value. 
This point is 207, which results in an error of one time sample, or a percentage error equal 

to }^5(,100% = 0.39% of Thop_min or approximately 20 ns. 





Figure 5.6: Signal in noise at SNR- 10 dB, ( a ) the ROC curve, ( b ) Pp vs. TfUreskoUi^ and ( c ) Pfa 
Ttkrtsk0id* Threshold, Tgureshou^ chosen as 30 times the variance of a detection vector generated in 
a *iio hop*^ frame. 
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D. RESULTS 


1. Detection 

The detection results given in Table 5.1 show the probability of detection, Pd, the 
probability of false alarm, Pfa, and the percentage of errors in classification for the 
selected threshold, Tthreshou, at each of the six SNRs considered. Note the entries under 
the column labeled ‘T’ represents a multiple of the variance of the detection vector 
generated from a “no hop” frame for each respective SNR level. The column labeled “% 
Error” shows the percentage of misclassifications (i.e., the percentage of false alarms plus 
misses). Note, also, that the a low probability of false alarm was sacrificed for higher 
probabilities of detection for reasons discussed earlier. For example, the entries in the 
row for SNR = 3 dB show that if a 11.4 % misclassification rate and a Pfa=0.1961 can be 
tolerated, then we can expect to detect 89.53% of the hops in a given frequency hopping 
signal. 


SNR 

k 

Pd 

Pf. 

% Error 

15 dB 

140 

1 

0 

0.0 

10 dB 

30 

0.9866 

0.0196 

1.4 

6dB 

15 

0.9844 

0.1569 

3.0 

3dB 

11 

0.8953 

0.1961 

11.4 

OdB 

1 

0.8129 

0.3529 

20.4 

-3dB 

3 

0.6927 

0.3333 

31.0 


Table 5.1: Detection statistics of 500 experiments 
applying the detection and estimation algorithm of 
this chapter. 
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2. Estimation 


Simulation results described in Section V-B are given in Table 5.2. The column 
labeled “Avg. Error” gives the average error obtained at each SNR level. For example, at 
the SNR level of 3 dB, out of all the hops which were detected, the average distance from 
the true hopping time was 10.48 sample points. This value equates to 
4.1 % of the minimum hopping time 2^^ • Columns with numeric headings indicate 

the hop detection probability within a given percentage of Thop mn- For example, at the 
SNR level of 3 dB, the column labeled “1%” indicates that 36% of all detected hops were 
located within 1% of Thop_mn or within 2 points of the true hop time, Thop- Similarly, 72% 
of all detected hops were located within 5% of Thop_min or within 12 points of the true hop 
time, Thop- Figure 5.7 plots the distribution of the hopping time detections for all SNR 
levels considered. Note that as the SNR decreases the distribution spreads out indicating 
less and less accuracy in the estimation, which is to be expected. 


SNR 

(dB) 

1% 

5% 

10% 

15% 

20% 

30% 

40% 

50% 

75% 

100% 

Avg. Error 
(# of samples) 

15 

0.790 

0.984 

0.992 

0.996 

0.996 

0.998 

1.00 

1.00 

1.00 

1.00 

2.22 

10 

0.726 

0.964 

0.974 

0.978 

0.982 

0.986 

0.986 

0.986 

0.986 

0.986 

2.70 

6 

0.558 

0.888 

0.926 

0.940 

0.950 

0.960 

0.968 

0.970 

0.970 

0.970 

5.46 

3 

0.360 

0.720 

0,758 

0.794 

0.828 

0.862 

0.874 

0.882 

0.886 

0.886 

10.48 

0 

0.II6 

0.302 

0.418 

0.510 

0.572 

0.684 

0.752 

0.768 

0.796 

0.796 

28.48 

-3 

0.090 

0.174 

0.276 

0.382 

0.456 

0.568 

0.614 

0.646 

0.686 

0.690 

30.99 


Table 5.2: Estimation statistics for the 500 experiments at each SNR level using the 
detection and estimation algorithm described in this chapter. Columns with numeric 
headings, show the probability that estimated hops are found within a given distance, 
expressed in percentage of Tkopjmm-> from true hopping times. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 

A. CONCLUSIONS 

This thesis considered the application of correlation functions and wavelet 
analysis to the detection and estimation of frequency hopping signals in additive white 
Gaussian noise. First, we described frequency hopping signals and briefly explained how 
they are used in spread spectrum communications techniques. Second, we introduced 
wavelet transforms and showed how they can be used to detect discontinuities in 
functions and their derivatives. Third, we introduced preprocessing techniques designed 
to improve the robustness of the detection and estimation scheme in noisy environments. 
Finally, the detection and estimation algorithm was presented and simulation results 
shown. 

The algorithm developed has only two restrictive assumptions: 

1. a minimum hopping time; 

2. a minimum frequency differential. 

Thus, it can find applications where the minimum hopping time is not held constant; i.e., 
in time hopping signals and hybrid techniques involving either/re^we/icv hopping or time 
hopping. Results show that the algorithm developed produces acceptable results for SNR 
levels as low as 3 dB. At this SNR level the percentage of misclassifications is 1 ].4% 
and the average error is 4.1% of the minimum hopping time 
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B. RECOMMENDATIONS 


Although the algorithm developed in this work succeeded in meeting its goals, it 
is believed that improvements could be obtained by altering the approach in two specific 
areas. First, the temporal correlation function produces a two-dimensional “image.” 
However, the algorithm developed here chose to perform one-dimensional wavelet 
analysis for each individual value of lag, t. Thus, a possible extension involves 
considering the problem as an image processing or a pattern recognition problem, due to 
the specific triangular pattern produced by calculating the TCF of frequency hopping 
signals. As a result, applying two-dimensional edge detection schemes, such as wavelet- 
based techniques, are expected to improve the robustness of the detection and estimation 
algorithm. In addition, one could use neural networks to recognize the triangular patterns 
of the TCF calculated for frequency hopping signals. Both extensions are presently under 
consideration. 
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APPENDIX A. MATLAB SOURCE CODE 
The Matlab source code for the detection and estimation algorithm developed in 
chapter five is provided in this appendix. Functions referenced, but not provided here are 
either a part of basic Matlab or a part of the Wavelet Toolbox [16], 


A. SIMULATION PROGRAM 


% This MATLAB code runs simulations of detection and 
9^ estimation of frequency hopping signals. The simulations 
have been designed to conduct 500 experiments at each of 
7 noise levels as listed in the SNR_db vector below, 

^ Uses DETECT2V, ERRLOC2 

^ Filename: d_esimul.m 

^ Capt Howard Overdyk, last revised 970902 

****•*#*«**«««4c«*4c*4c««9|c*3)c*»e************************************** 

9^ Create required signals: 

clear 


load seed, % Using set seed only to ensure results are 

rand('sced\seed); 9^ rcproduceable. 

dcllaJ=IOOO. 

T_min=256: 

Cmin=IO^; 
f_max=24*IO^; 
fs=50*IO^; 

Nf=(f_max-Lmin>/dclia_f; 

Nc-MK). 


% 1 KHz 

% T_min = 256 pts = 5.12 us 
% 1 MHz 
% 24 MHz 

9^ Sampling frequency = 50 MHz 
9^ Number of different random freqs fk 
9^ Number of experiments 


^ Prixlucc random hop times, Th. 

Th=round(206*rand( l.Ne)); 9^ 206 to allow for removing 25 pis from 

9^ each end. 0 => no hop 
9c In practice, overlapping technique 
9r would be employed to account for 
9c removing edges. 

for i2=0:49. 


Th( l+i2*10)=0; 9c ensure 10^ are no hops to fully test 

9c detection portion of algorithm 
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end 


rand('seed’,seed); 

% Produce random frequencies, Fh, where 1st row is hop from freq. 

% and the 2nd row is the hop to freq. 

Fh=(round(Nf*rand(2,Ne)).*delta_f)+f_min; 

hops=zeros( 1 ,Ne); % to keep track of if there is a hop and where, 

% zero (0) => no hop or else it will be hop pt. 


x=[]; % Initialization purposes 

% Use random hop times, Th, and random freqs, Fh, to produce signals 
% for 500, Ne, experiments. 

for iI=l:Nc, 

ifTh(iI) = 0, 


n=l:T_min; 

x=[x cos(2*pi*Fh(l,il)/fs*n)]; 
clseifFh(l,il) = Fh(2,il), 


n=l:T_min; 

x=[x cos(2*pi*Fh(l,il)/fs*n)]; 


else 

hops(il)=Th(il)+25; % RN btwn 26 and 231 

n 1 = 1 :hops(i 1); % for first freq 

n2= I :T_min-hops(i!); % for second freq 

x=IX cos(2*pi*I^( 1 .i I )/fs*n 1) cos(2*pi*Fh(2,i 1 )/fs*n2)]; 

end ^ if statement 

end 91 loop 

91 

91 Create Noise and embed signal in it and make analytic. 

SNR_db=[ 15 10 6 3 0 -3); 9c SNR in dB 

Thresh=:1140 30 15 11 I 3); % Thresholds pre-determined from ROC 

% function 

sigma=(sqrt(I0.^(-SNR_db( 1 )/10)./2)... 
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sqrt(10.''(-SNR_db(2)/10)./2)... 
sqrt(10.''(-SNR_db(3)/10)./2)... 
sqrt(10.''(-SNR_db(4)/10)./2)... 
sqrt(10.''(-SNR_db(5)/10)./2)... 
sqrt(10.^(-SNR_db(6)/10)./2)]; 


randn('seed',seed); 
w 1 =randn(6,length(x)); 

for i=l:6, 

w(i,:)=sigma(i).*wl(i,:); 

xn(i,:)=hilbert(x+w(i,:)); 


end 

xn(7,:)=hilbert(x); % No noise case 




% 


% Perform Experiments 


fori 1 = 1:7, 

for i2=l:Ne, 


ld( i I ,i2),d(i 1 +7,i2)]=delect2v(xn(i 1 ,(i2-1 )*T_min+1 :i2*T_min)); 
end 


Thrcshld=Thresh(il)*d(il+7J); % First experiment is known 

% to be a no hop case 

indl=d(il+7,:)>Threshld; % Create an indexing vector 

Which equals one only if the 
^ variance of the detection 
9^ vector exceeds the pre- 
^ determined threshold. 

da=d(il,:) *indl; % All detections which do not 

^ exceed the threshold are set 
% to a miss 


»«i**«***•***»«****♦**««**«««*«*««4t4t**«*««*«**»«««♦ 

91 

91 Analysis of Results call crrloc2 


|dclla_th(i 1 ,:).Pd(i 1 ),Pfa(i 1 ),pcnl_crr(i 1 ),avg_err(i 1 ))=errloc2(da,hops,Fh); 
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B. DETECTION FUNCTION 


function [t,v]=detect2v(x) 

% [T,V]=DETECT2V(X) given signal X, detects the hop time of 256 
% samples of a signal X and returns the hop point in T and 
% the var(WT summed coeffs) in V. 

% 

% Uses TCF3, TRDET3U3, SUMLAGS3 
% Filename: detect2v.m 


% Capt Howard Overdyk, last revised 970813 




min_lag=0; 

max_lag=length(x)/2; 

Rxt=tcD(x.min_lag,maxJag); 

c=trdet3u3(Rxt,min_lag,maxJag); 

(cal,ca2,ca]=sumlags3(c); 

cs=(0 ca( 1 :length(ca)-1 )]+[ca(2:length(ca)) 0]; 


v=var(cs); 

(in,tl=max(cs); 


C. TEMPORAL CORRELATION FUNCTION 

function Rx=tcf3(x» arg2,arg3) 

TCF3{XARG2,ARG3) returns the temporal correlation funciionfTCF) 
9^ of X. If arg2 is present the tcf will only be calculated 

from min_lag=arg2 out to max_lag=arg3, otherwise the whole 
TCF will be computed for positive values of lag, tau. 

Filename: lcf3.m 

9^ Capt Howard Overdyk, last revised 970512 








if nargin>l 

minjag=arg2; 

maxjag=arg3; 

else 

minjag=0; 

inaxJag=length(x)/2; 

end 

Rx=zeros(length(x)+1 ,length(x)/2+1); 
for lag=minjag:maxjag; 
for i=lag+l:length(x)-lag, 

Rx(i,lag+l)=x(i+lag)*conj(x(i-'lag)); 

end 

end 

% Compute the TCP for negative values of lag, tau. 

%Rx(:,l:length(x)/2)=conj(fliplr(Rx(:,length(x)/2+... 
% 2:Iength(x)+l))); 


D. WAVELET TRANSFORM OF TCF FUNCTION 

function ca=trdct3u3{R,minJag,maxJag) 

% TRDET3U3(R,MIN_LAG,MAX„LAG) performs wavelet transforms, using 
% Haar wavelet, on all lags from MIN_LAG to MAX_LAG on phase 
% of temporal correlation fen, R, and then sums the 
9^ coefficients of the first two scales of each transform. 

9f Returns the matrix of summed coefficients CA. 

9f Uses UNWRAP. ANGLE. DIFF, MEORLT, WLOOKNP 

9f Filename trdct3u3.m 

Capt Howard Overdyk. last revised 970420 

Extract unwrapped angle of TCF. 

Rxa-un wrapt angle( R)); 
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% Investigate different lags of TCF of complex signal to 
% determine hop time 

level=4; 

ca=zeros( 129,length (R(:, 1))-1); 
for lag=min Jag:max_lag, 

% median filter Rxa with length 5, then differentiate 
% and finally, median filter with length 25 before 
% performing the wavelet transforms. 

Rxtmf=medfilt(diff(medfilt(Rxa(:,lag+l).',lag,5)),lag,25); 
c=wtnp(Rxtmf,'dbr,level);% Haar == dbl 

ca(lag+l,:)=sum(abs(c(l:2,:))); % sum the 1st two scales 


end 


E. MEDIAN FILTER FUNCTION 

function xmf=medfilt(x,lag,f_size) 

% MEDFILT(X) given data X, which LAG the data came from, and 
% filler size, F_SIZE, this function returns the median 
% filtered data in the vector XMF. 

% LAG is used to determine how large of area outside of 
% the TCF is zeros and so, the function can set these 

% values to the median of the first three values within 

% the TCF. 

^ Uses MEDIAN 

% Filename: medfilt.m 

% 

^ Capt Howard Overdyk, last revised 970421 

91 

if nargin < 3 

Lsize=r3; 


end 

for i2=l:lag. 
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x(i2)=median(x(lag+l :lag+3)); 

x(length(x)-i2+l)=median(x(length(x)-lag-2:length(x)-lag)); 

end 

for i= 1 :ceil(f__size/2), 

xmf(i)=median(x( 1 :f_size)); 
end 

for i=ceil(f_size/2)+l :length(x)-ceil(f_size/2), 

xmf(i)=median(x(i-ceil(f_size/2)+i:i+floor(f_size/2))); 

end 

for i=:length(x)*ceil(f_size/2)+l :length(x), 

xmf(i)=median(x(length(x)-f„size+l:length(x))); 

end 


F. FOUR LEVEL WAVELET TRANSFORM FUNCTION 

function m=wtnp(x,wavfamn,level) 

% WTNP(X,WAVFAMNLEVEL) given signal, X, 

% WAVFAMN, the wavelet name (e.g. ’db2’), and 
% the dyadic LEVEL for decomposition, the function 

% returns the wavelet transform coefficients in 

% matrix, M. 

% 

% Level <= 4 

^ Uses COEFMAT4, WAVEDEC 

Filename: wtnp.m 

% Capt Howard F. Overdyk, last revised 2/7/97 

^mm*mm*********mmmmm************************************* 


lc.l)=wavedcc(x,icvel,wavfamn); 


m=cocfmai4(c,l,level,wavfamn); 
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G. WAVELET TRANSFORM COEFFICIENT MATRIX FUNCTION 

function data=coefmat4(c,I,Ievel,arg3,arg4) 

% COEFMAT4(C,L,LEVEL,arg3,arg4) returns an level x Ns matrix of 
% the wavelet coefficients (details only) Unnormalized 
% coefficients using the DWT tiling and 
% using either a specific wavelet (‘wname’ see WFILTERS) 

% or specific wavelet decomposition filters. 

9c Ns is the length of the signal X 
% 

9c The structure is organized as: 

9( C - [app. coef(N)ldet. coef(N)l... Idet. coef.(l)] 

9c L( 1) = length of app. coef.(N) 

L(i) = length of det. coef.(N“i+2) for i = 2,...,N+1 
L(N+2) = )ength(X). 

9c 

9 Level <= 4 

9 

9 Uses DETCOEF, WFILTERS, WKEEP 

9 Filename: coefmat4.m 

9 

9 Howard F. Overdyk last modified 970807 

^ 0»m0m*»*»***mm)$m**********3t^************************************** 

9 

9 Check arguments 
if nargin=:=4 

|LoF_D,HiF_D]=wfiliers(arg3,'d'); 

else 

LoF_D = arg3; HiF_D = arg4; 
end 

9 Extract approximation and details 

N=l(lcngth(l)). 

9 Expand d I 

d I n=dclc«f(c J, I). 

dl=wkecp<dln,N/2^LT); 

il=:!:2:N.i2=2:2:N: 
init=zeros( LN); 
ddl(il>^l;ddl(i2)=^l; 

9f Expand d2 
if (level >= 2) 
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d2n=detcoef(c,l,2); 

d2=wkeep(d2n,N/2''2;i'); 

il=l:4:N; i2=2:4:N; i3=3:4:N; i4=4:4:N; 

dd2=init; 

dd2(i 1 )=d2;dd2(i2M2;dd2(i3)=d2;dd2(i4)=d2; 
end 

% Expand d3 

if (level >= 3) 

d3n=detcoef(c,l,3); 

d3=wkeep(d3n,N/2'^3,'l'); 

il=l:8:N; i2=2:8:N; i3=3:8:N: i4=4:8:N; 

i5=5:8:N; i6=6:8;N; i7=7:8:N; i8=8:8:N; 

dd3=init; 

dd3(i 1 )=d3;dd3(i2)=d3;dd3(i3)=d3;dd3(i4)=d3; 
dd3(i5)=d3;dd3(i6)=d3;dd3(i7)=d3;dd3(i8)=d3; 
end 

% Expand d4 

if (level >= 4) 

d4n=detcoef(c,l,4); 

d4=wkeep(d4n,N/2M,'r); 

il=l:16:N; i2=2:16:N; i3=3:16:N; i4=4:16;N; 

i5=5:16:N; i6=6:l6;N; i7=7:16;N; i8=8:16:N; 

i9=9:16:N;il0=10:16:N;ill=ll:16:N;il2=12:16:N; 

il3=13:16:N; il4=14:16:N; il5=15:16:N; il6=16:16:N; 

dd4=init; 

dd4(i 1 )=d4;dd4(i2)=d4;dd4(i3)=d4;dd4(i4)=d4; 
dd4(i5)=d4;dd4(i6)=d4;dd4(i7)=d4;dd4(i8)=d4; 
dd4(i9)=d4;dd4(ilO)=d4;dd4(ill)=d4;dd4(il2)=d4; 
dd4(i 13)=d4;dd4(i I4)=d4;dd4(i 15)=d4;dd4(il6)=d4; 
end 


dala=Iddl:dd2;dd3:dd4]; 

H. 45“ MATRIX SUMMATION FUNCTION 

function |cal,ca2, ca]=sumlags3(c) 

9i SUMLAGS3(C) performs a 45 degree summation over all columns of the 
9^ wavelet transform coefficients, C, in both directions 
9^ of 45 deg and returns a row vector CA with the result. 

9^ CAl contains results of summing only in 45 deg direction. 

9r CA2 contains results of summing only in 135 deg direction. 

Edges arc clipped to avoid edge effects using ERLTZ. 

9c Uses results of EDGEFILT2 

9c Filename: sumlags3,m 

9c 

9c Capl Howard Ovcrdyk, last revised 970707 
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Steele********************************************************* 

% 


load efilt2; 

[m,n]=size(c); 

c=efilt2.*c; 

cal-zeros(l,n); 

ca2=zeros(l,n); 


% see code of edgefilt2.m 
% Remove edges of WT coeffs. 


for il = l:n-m+l, % covers area where we have all rows to 

% sum over at 45 degrees 


for i2= l:m, 


cal (i 1 )=cal (i 1 )+c(i2,i 1-1 -f-i2); 


end 

end 

cal=cal./m; % Normalize summation 


ci=l; 


for il = n-m+2:n, % covers end where we aren't summing over 

% all rows at 45 degrees. 


for i2 = l:m-cU 

ca 1 (i 1 )=ca 1 (i 1 )+c(i2,i 1 -1+i2)./(m-c 1); 


end 

cl=cl+l; 


end 


for i I = n:-1 :m, % covers area where we have all rows at 135. 


for i2 = m:-l :1, 

ca2(i I )=ca2(i 1 )-K:(i2d 142+1); 


end 

end 

ca2=ca27m; % Normalize summation 


cl = l; 


84 



for i 1 = m-11:1, % covers area where we don't have all rows 

% at 135 degrees. 


for 12 = 

ca2(i 1 )=ca2(i 1 )+c(i2,i 1 -12+1 ),/(m-c 1); 


end 

cl=cl + l; 


end 

ca=cal+ca2; 

I. EDGE FILTER CREATION PROGRAM 

% EGDEFILT creates a matrix which will filter off the edge effects 
% of a wavelet transform of a TCP of size 129 x 256. 

91 Filename: edgefilt2.m 

91 Capi Howard Overdyk, last revised 970707 

^ •*«**«* ^ »****# )|e:|c * ifesie # 9)c 9f: ** 9ic)ic:f: 9|c 9(; 4: *** 9^ ********* >|( Nc * 4: ^ ***)»: ****** 4c 

clear 

cuioff=25; 

m=l29; 

n=256; 

mis 128-cutoff; 

. fl=zcros(mKml); 

for 

f I =f 1 +diag(ones( 1 j I ),(i 1 -m 1)); 


end 

f2=|niplr(fl)fl|; 

f2=nipud(f2); 

f=|r2; zcros(26.2*ml)l; 

crili2=Izcros(m.25) f zcros(m.25)l; 

save /homc/dsp()b/overdykymallab/ihcsis/fq_wk7/erill2.mal efill2; 

J. DETECTION AND ESTIMATION ANALYSIS FUNCTION 

function (dclta_th,Pd.Pfa,pcnt_crr.avg_crT)=eiTloc2(d,hops) 
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% [DELTA^TH,Pd,Pfa,PCNT3RR,AVG_ERR]=ERRL0C2(D,H0PS) will 
% given: D - a 1x500 vector containing the location of the 
% detected hops for each experiment. 

% HOPS - a 1x500 vector containing the actual location 

% of the hops for each experiment. 

% 

% returns 
% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% side effects: 

% Produces two plots of the distribution of the location 

% of the detections relative to the actual hopping time 

% once in line graph format and then in bar graph format. 

% Filename: errloc2.m 

% 

% Capt Howard Overdyk, last revised 970805 

% 

% Analysis of Results 

Ne=500; 

Tol_hops2=Nc; 

Ntcf=256; 

N=Nc*Nicf; 
max_dih=Nicf» 
dih_ ind=0:max_dih. 
dclia_lh=zcros( I .max_dih+1); 
dclia_in=zcros( I ,max.dih+1); 
loc„ind=0:Ntcf; 
loc=zcros( 1,257): 

LoLni=zcros{2.Nc); 9^ first two rows correspond to freq either 

side of missed hops on experiment # = col # 

9^ Initialize output variables 


DELTA_TH- a 1x257 vector containing the distribution 
of correctly detected hops versus the 
distance from the actual hop time in # 
of samples. 

PD ” a scalar value of the probability of detection 

PFA - a scalar value of the probability of a 

false alarm. 

PCNT_ERR- a scalar value of the percentage of mis- 
classifications. 

AVG_ERR - a scalar value indicating, of the actual 
detections, what the average distance, in 
# of samples, from the true hop time was. 






dets=Pd; 

fa=:Pd; 

miss=Pd; 

nn=Pd; 

i 1 = 1; % remains constant 

for i2=l :Ne, % the number of experiments 

acth=hops(i2); 

exph=d(i 1 ,i2); % hops detected by experiment 

delta_h-abs(acth>exph); 

if acth ~= 0, % hop exists 

if exph 0, % found hop 

delta_th(i 14elta_h+1 )=delta_th(i 1 ,delta_h+l)+1; 

delta_m(i 1 ,i2)=delta_h; 

dets(il)=dets(il)+l; 

else % missed detection 

loc(i 1 ,acth)=loc(i 1 ,acth) + 1; 
miss(il) = miss(il)+ 1; 
end % found hop if statement 

else % no hop exists case 

if exph ~= 0» % false alarm case 

fa(il) = fa(il)+ 1; 

else % no hop no detect case 

dclla_th(i 1 .delta_h+1 )=della_th(i 1 ,delta_h+1)+1; 
della_m(i 1 ,i2)=delta_h; 

9r della^h always equals zero (0) here. 
nn(il) = nndIHl; 

end ^ false alarm if statement 

end 9^ hop exists if statement 

end 9^ Inner for loop i2 
9^end 9^ Outer for loop il 
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figure 
orient tall 


[i3,i4,v]=find(hops--=0); 

% Print relavent statistics to MATLAB session 

Toi_hops=sum(v) 

dels 

miss 

fa 

nn 

Pd=dets/Tot_hops 

Pfa=fa/(Ne-Tol_hops) 

subplol(22I) 

plot(dlh_ind,delta_th( 1, :)/Tot_hops2,’m-’); 
xlabel('dellaTh'),ylabel('Detections') 
titieCDislribution Hopping Time Detection’) 

^ prim out in tabular form the percentage of detections 
^ within 1,5, 10, 15, 20, 30, 40,50,75, and 100 percent 
^ of the minimum hop time. 

p_tolhop=sum(della_th.’)/Tot_hops2 
p_iot l=sum(dclia_th(:, 1:3).’)/Tot_hops2 
p jot5=sum(della_th(:, 1:13).')/Tot_hops2 
p_tot 10=sum(delta_th(:, 1:26).’)/Tot_hops2 
p_tol 15=sum(dclla_lh(:, 1:39).’)/Tot_hops2 
p_tol20=sum(delu_th{1:52).’)/Tot_hops2 
p_iot30=sum(dclta_th(:,l:77).')/Tot_hops2 
p_tot40=sum(dclia_lh(:, 1:103).')/Tot_hops2 
pjot50=sum(delia_lh(:,l: 129).')/Toi_hops2 
p_tot75=sum(dclla_lh(:, 1:193).*)/Tol_hops2 
p^iol 100=sum(dclta_th.')/Tot_hops2 


subploi(222) 

bar(0:max_dlh.dclia_lh( 1 ,:)/Tot_hops2); 
axis(10.max_dth,-inf,infl) 
titleCDistribution Hopping Time E)cleciion’) 

a vg_cm=mcan(della_m) 
pcnt_crr=(fa + miss)/Nc 
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K. RECEIVER OPERATING CURVES GENERATION FUNCTION 


function [Pd,Pfa]=roc(d,hops, inc_size) 

% [Pd,Pfa]=ROC(D,HOPS,INC_SIZE) given a vector, D, of detected hops 
% a vector, HOPS, of actual hops, and INC_SIZE which is the 

% amount to increase the threshold by at each increment, will 

% return a probability of detection vector, PD, and a 
% probability of false alarm vector, PEA. 


% 

% 


Filename: roc.m 

Capt Howard Overdyk, last revised 970913 


% 

% Analysis of Results 


Ne=500; 

Tot_hops2=Ne; 

Ntcf=256; 

N=Ne*Ntcf; 
max_dth=Ntcf; 
dth_ind=0:max„dth; 
delta_th=zeros( 1 ,max_dth+1); 

^ Initialize output variables 

Pd=zeros( 1,100); 

Pfa=Pd; 

dcls=Pd; 

fa=Pd; 

miss=Pd; 

nn=Pd; 


Th=d(2,l) 

|i3j4,v)=find(hops-=0); 

Toi_hops=sum(v) 


foril = l:IOO; 

TTh=inc_size*i 1 *Th; 

indl=d(2.:)>TTh; 

dl=d(l.:).*indl; 

for i2= 1 :Ne, % the number of experiments 

acth=hops(i2); 
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exph=dl(i2); 


% hops detected by experiment 


delta_h=abs(acth-exph); 

if acth -= 0, % hop exists 

if exph 0, % found hop 

dets(i 1 )=dets(i 1)+1; 

else % missed detection 

miss(il) = miss(il) + 1; 
end % found hop if statement 

else % no hop exists case 

if exph 0, % false alarm case 

fa(il) = fa(il)+ 1; 

else % no hop no detect case 

% delta_h always equals zero (0) here. 

nn(il) = nn(il)+l; 

end % false alarm if statement 

end % hop exists if statement 


end % Inner for loop i2 

end % Outer for loop i 1 

Pd=dcts/Tot_hops; 
Pfa=fa/(Ne-Tot_hops); 


figure 

subploi(22l) 

ploi(Pfa.Pd) 

xlabel(*Pfa*).ylabcl(Td‘) 

lillcCROC) 

grid on 

k=inc_sizc*[ 1:100]; 

subploi(222) 

plol(k.Pd) 

litlcCPd vs Threshold multiple*) 
grid on 



subplot(224) 

plot(k,Pfa) 

title(Tfa vs Threshold multiple') 
grid on 







APPENDIX B. THRESHOLD DETERMINATION 


Figures B.l to B.5 plot the curves used to select detection threshold values for 
SNR levels of 15 dB, 6 dB, 3 dB, 0 dB, and -3 dB, respectively. The threshold, Tthreshoid, 
is given by: 

T,hresho!d = ^ ’ ^^^[[detection vector^j^p{t)) , (B.l) 

where k is the multiple which needs to be determined. The horizontal axis in plots ( b ) 
and ( c ) of Figures B.l to B.5 represents the parameter k. Therefore, the selection of k is 
based on choosing an acceptable Pd level or an acceptable PfA level. The parameter, k, is 
chosen so that it leads to a high probability of detection, Pd, and an acceptable Pfa level. 
This selection procedure is explained in further detail in Section V.C.7. 

For example, let us consider the 6 dB case plotted in Figure B.2. Let us assume 
that we need Pd = 0.985. Figure B.2b shows that the corresponding required value of k 
equals 15. This value of k results in a Pfa level equal to 0.16. Note, that these values of 
Pd and Pfa are located at a point very near the “elbow” of the RCXT curve shown in 
Figure B.2a. 
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ROC Pd vs Threshold mu K pie 



Figure B.l: Noisy signal with SNR = 15 dB, (a ) ROC curve, ( b ) Pd vs. and ( c ) Pfa vs. 

TutmhM- Threshold, T^miuu^ chosen as 140 times the variance of the detection vector generated in a 
*'no hop” frame. 
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Figure B.2: Noisy signal with SNR s 6 dB, ( a ) ROC curve, ( b ) Pp vs. Tu^tskmU^ ( c ) Pfa vs. 
Ttkntk 0 U* Threshold, T^^ntkou^ chosen as 15 times the variance of the detection vector generated in a 
**00 hop” frame. 
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Figure B.3: Noisy signal with SNR = 3 dB» (a ) ROC curve, ( b ) Pp vs, Tthmhou^ and ( c ) Pfa vs. 
TthrrtMd^ Threshold, Ttkreskou^ chosen as 11 times the variance of the detection vector generated in a 
hop” frame. 
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Figurr B.4: Noisy signal with SNR = 0 dB, ( a ) ROC curve, ( b ) Pp vs. and ( c ) Pfa vs. 

^rtrrrrrfi^ Threshold, chosen as 1 times the variance of the detection vector generated in a ^Tio 

hop** frame. 
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Figure B.5: Noisy signal with SNR = -3 dB, ( a ) ROC curve, ( b) Pd vs. Toutthou^ and ( c ) Pfa vs. 
TArttkM- Threshold, T,m»U’> chosen as 3 times the variance of the detection vector generated in a “no 
bop” frame. 
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